In [1]:
# Setting up a model and a mesh for the MT forward problem
In [24]:
import SimPEG as simpeg
from SimPEG import MT
sys.path.append('/home/gudni/gitCodes/python/telluricpy')
import telluricpy
In [6]:
# Define the area of interest
bw, be = 556500, 558000
bs, bn = 7133500, 7133500
bb, bt = 0,480
In [14]:
# Build the mesh
# Design the tensors
hSize,vSize = 50., 12.5
nrCcore = [10, 8, 6, 4, 2, 2, 2, 2, 2]
hPad = simpeg.Utils.meshTensor([(hSize,9,1.5)])
hx = np.concatenate((hPad[::-1],np.ones(((be-bw)/hSize,))*hSize,hPad))
hy = np.concatenate((hPad[::-1],np.ones(((bn-bs)/hSize,))*hSize,hPad))
airPad = simpeg.Utils.meshTensor([(vSize,13,1.5)])
vCore = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcore,(simpeg.Utils.meshTensor([(vSize,1),(vSize,8,1.3)])))])[::-1]
botPad = simpeg.Utils.meshTensor([(vCore[0],8,-1.5)])
hz = np.concatenate((botPad,vCore,airPad))
# Calculate the x0 point
x0 = np.array([bw-np.sum(hPad),bs-np.sum(hPad),bt-np.sum(vCore)-np.sum(botPad)])
# Make the mesh
meshFor = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
In [17]:
print np.sum(vCore)
print meshFor.nC
print meshFor
In [35]:
# Save the mesh
meshFor.writeVTK('nsmesh.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh.vtr')
In [ ]:
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)
In [39]:
telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)
In [40]:
# Get active indieces
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')
In [83]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
'XVK':3e-2,
'PK1':5e-2,
'PK2':1e-2,
'PK3':1e-2,
'HK1':1e-3,
'VK':5e-3}
In [75]:
# Loop through
extP = '../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')
In [89]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','XVK','PK1','PK2','PK3','HK1','VK']:
sigma[geoStructIndDict[key]] = geoStructFileDict[key]
In [90]:
# Save the model
meshFor.writeVTK('nsmesh_0.vtr',{'S/m':sigma})
In [91]:
# Set up the forward modeling
freq = np.logspace(5,0,26)
np.save('MTfrequencies',freq)
In [58]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw,be,bs,bn,bb,bt)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=True)
In [92]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+25.,be,50),np.arange(bs+25.,bn,50))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,modTopoVTP)
np.save('MTlocations',locArr)
In [60]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py
In [95]:
%matplotlib qt
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf
In [99]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_0.npy')
mtData
Out[99]:
In [110]:
iPf.MTinteractiveMap([mtData])
Out[110]:
In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions,
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.
Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py
In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')
In [109]:
np.unique(drecAll['freq'])[10::]
Out[109]:
In [ ]: